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Abstract 

We study a network of 500 coupled modified van der Pol oscillators. The value 
of a parameter associated with each oscillator is drawn from a normal distribu- 
tion, giving a heterogeneous network. For strong enough coupling the oscillators 
all have the same period, and we consider periodic forcing of the network when 
it is in this state. By exploiting the correlations that quickly develop between the 
state of an oscillator and the value of its parameter we obtain an approximate 
low-dimensional description of the system in terms of the first few coefficients in a 
polynomial chaos expansion. Standard bifurcation analysis can then be performed 
on this low-dimensional system, and the results obtained from this predict very 
well the behaviour of the high-dimensional system for any set of realisations of the 
random parameter. Situations in which the method begins to fails are also discussed. 
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1 Introduction 



Synchronisation is a common phenomenon in biology and elsewhere [T|15!l25j . 
It is often studied by investigating the conditions under which oscillators in a 
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particular network will synchronise [3][2Tj . Periodic forcing of systems is also 
ubiquitous [15], and so it is natural to study the entrainment of a network of 
coupled oscillators by a periodic forcing function. Many authors have studied 
small networks of two or three non-identical oscillators [3] , and larger networks 
of oscillators that may have some symmetry [T6f2T] or a particular form of 
coupling [29J . The continuum limit in which there exists an infinite number of 
oscillators has also been studied in detail and many results are known for this 
case [2fHl33] . However, it is known that finite networks can show behaviour 
that does not occur in the continuum limit [UTTT] . In many situations, finite 
networks are the most realistic way to model a physical system [5llll|l30| . 
Results for large, finite networks will thus help bridge the gap between small 
network dynamics (for which bifurcation analysis is straightforward) and those 
for an infinite number of oscillators (where statistical physics provides the 
appropriate tools). 

In this paper we consider a large but finite heterogeneous network of coupled 
oscillators, which are collectively periodically forced. However, we do not anal- 
yse the system exactly; instead we analyse a low-dimensional description of it. 
This is the "equation- free" approach developed by Kevrekdis et al. [20] . The 
results here extend those of Moon et al. [26i 28] : we consider two- variable oscil- 
lators, capable of undergoing Hopf bifurcations; we consider periodic forcing of 
the network, and we perform bifurcation analysis on the system to understand 
how the behaviour of the system changes as parameters are varied. 

The system we study is 



dt ~ m Xl 



iV 3=1 



^j- = —Xi + A sin (cut) (2) 

for i = 1,...,N, where N is the number of oscillators in the network. For 
most of this paper we set N = 500. The oscillators are van der Pol oscil- 
lators [IT] with an extra term (xf/2) which breaks the internal symmetry 
[(x, y) — > (—x, —y)\ of each individual oscillator. These oscillators were chosen 
as being "typical" in the sense of not having any particular properties. For 
(3 = 0, an uncoupled oscillator (e = 0) undergoes a supercritical Hopf bifur- 
cation as <p increases through zero, with angular frequency 1. The angular 
frequency for an isolated oscillator as a function of is shown in Fig. [TJ 

The Hi are taken from a normal distribution with mean and standard devia- 
tion 1. (As discussed below, the methodology can be used with other distribu- 
tions.) If (3 the network is heterogeneous, and each oscillator, if uncoupled, 
would have a different angular frequency determined by the value of + (3^. 
When A = 0, for f3 small enough, <fr of moderate size and e large enough, 
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the oscillators synchronise in the sense of having the same period. Note that 
oscillators % and j cannot synchronise in the sense of Xi(t) = Xj(t) for all t 
unless fa = fa. In this synchronised state the attractor of the system is a 
periodic orbit, which could be parametrised by a periodic variable, say 6{t). 
The variables xi, . . . xn, yi, ■ ■ ■ , Un could each then be written as functions of 
9. This description would no longer be valid if one or more of the oscillators 
"unlocked" from the group. 

We want to study the system in this synchronised state, but do not want to 
keep track of all the 2N variables x\, . . .xn, yi, . . . , yjy. Instead, we describe 
the state of the system by a small number of variables. We cannot easily derive 
an equation that governs the dynamics of these variables, but by repeatedly 
mapping between the two levels of description of the system we can numer- 
ically evaluate the results of integrating these unavailable equations; we can 
also find their collectively periodic states and their dependence on parameters, 
without ever obtaining the reduced equations in closed form. 

If the system is in this synchronised state and we increase A from zero, it will 
become periodically driven and it may be possible for the oscillators to lock 
with the driving frequency [15]. The latter part of this paper will consider 
this phenomenon in detail, but we first discuss the particular low-dimensional 
description of the forced system ([I])-© used here, and how it can be used in 
projective integration to speed up direct simulation of this system. 



2 A low-dimensional description 

The main idea behind the low-dimensional description used here depends on 
correlations that rapidly develop between x iy yi and the value of fa in the 
parameter regime where synchronization eventually prevails. This is demon- 
strated in Fig. [2] where we plot the Xi and yi as functions of \ii for a particular 
realisation of the fa at three different times. We see that after just two peri- 
ods of the forcing strong correlations develop between the state of an oscillator 
and its fa value. We will see that these correlations occur whether or not the 
network is synchronised with the forcing. 

These correlations allow us to expand the x and y in certain classes of polyno- 
mials of \i [26] ; Hermite polynomials are appropriate for a normal distribution 
of fa while different /i distributions correspond to different polynomial types 
(the so-called Generalised Polynomial Chaos, or GPC [M]). We write 

q 

x(t,fa) = Y,aj(t)H j (fa) (3) 

j=0 
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(4) 



3=0 



where Hj is the jth Hermite polynomial [H (x) = l,Hi(x) = 2x,H 2 (x) = 
Ax 2 — 2, . . .]. This expansion is known as a polynomial chaos expansion [26J, and 
the aj and fcj are the polynomial chaos coefficients. For a specific realisation of 
the fa, we have Xi(t) = x(t, fa) and similarly for the y%{t). Our low-dimensional 
description then involves the 2(q + 1) coefficients a , . . . a q , b , . . . , b q . This de- 
scription is approximate, and the approximation becomes better as q is in- 
creased. Given Xi and z/t for a particular set of //j, the are found by min- 
imising the quantity 



N 

E 

i=l 



(5) 



3=0 



and the fej are found by minimising 

2 

W -Yl h 3 H 3^i 

3=0 



N 

E 



(6) 



This is easily done in Matlab using the "backslash" operator to solve an overde- 
termined linear system. The operator from the Xi and yi to the dj and 6j is 
referred to as the "restriction" operator. Similarly we construct a "lifting" 
operator: given the di(t) and bi(t) and a particular realisation of the \i% we 
have 



i 

X i (t) = J2 a j( t ) H j(l J 'i) ( 7 ) 
3=0 

y l {t) = Y,b J {t)H J {H % ) (8) 

3=0 

Armed with these two operators we can now proceed to numerically solve the 
unavailable equation for the polynomial chaos coefficients. 



3 Coarse Projective Integration 



Coarse projective integration entails accelerating the simulation of a high- 
dimensional system by projecting forward in time using only the variables 
in a low- dimensional description of the system. This is accomplished by oc- 
casionally performing short bursts of full simulation of the high-dimensional 
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system in order to obtain the numerical information (estimation of the time- 
derivatives of the low dimensional description variables) required to perform 
accurate projections [2"Uj26f28| . We can use the low-dimensional description 
in the previous section for coarse projective integration as follows. For conve- 
nience, let the high-dimensional description be the variable 

X =[x 1 ,...,x N ,y 1 ,...,y N }eM 2N (9) 

and the low- dimensional polynomial chaos coefficient description be the vari- 
able 

Z= [a ,...,a q ,b ,...,b q ] Gl 2(9+1) . (10) 

Given X(0), integrate ([I])-© forward for Ni steps of size St. Calculate Z at 
some or all of the times t = 0, St, 2St, . . . , N\St using the restriction opera- 
tor. Use these values of Z to extrapolate the values of Z to a time N 2 St in 
the future, i.e. to time (Aq + N 2 )St. Lift from the value of Z{{N X + N 2 )St) 
to X((Nx + N 2 )St) as detailed above. Restart the integration of ([I])-© us- 
ing X((Ni + N 2 )St) as the initial condition and integrate for a further Aq 
time steps. Restrict to Z and repeat the procedure. If the cost of restrict- 
ing, extrapolating and lifting is small compared to the cost of integrating the 
system (CO)-® for time steps, this procedure may well be faster than in- 
tegrating ([I])-© directly. We expect this to be the case when the full system 
is characterized by a separation of time scales; the same principle underpins 
several analytical reduction techniques (e.g. centre manifolds, inertial mani- 
folds) but in our case the reduction is obtained on the fly, from the short full 
simulation bursts. 

We show results in Figs. [3] and H] for St = 0.005, Aq = 3. In this case the 
projective step involves fitting a cubic to the last Aq + 1 data points (the 
last Aq of which were obtained through direct integration of the full system) 
and then evaluating this cubic at a time N 2 St in the future. The top panel 
of Fig. [3] shows the speedup as a function of N 2 . The speedup is defined as 
the time taken to directly integrate ([I])-© over < t < 100 with time-step 
St divided by the time taken to integrate over < t < 100 using coarse 
projective integration, as described. A speedup greater than 1 (N 2 greater 
than approximately 10) means that projective integration is more efficient 
than direct integration (provided that the accuracy in the values of the reduced 
variables is satisfactory). 

Of course, as N 2 is increased the integration will start losing accuracy. The 
bottom panel of Fig.[3]shows the results of an integration when N 2 = 1 and the 
two curves shown (one for projective integration and one for direct integration) 
are indistinguishable. The top panel of Fig. H] shows the case when N 2 = 71. 
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In this case, coarse projective integration involves taking 3 steps of length St, 
giving the clusters of 4 points shown in the bottom panel of Fig. HJ and then 
projecting the a, and foj forward a time 7lSt, lifting these values to initialise 
the Xi and and continuing. The integration is clearly less accurate than that 
shown in Fig. [3], but the general behaviour is still qualitatively reproduced. 
In the spirit of Taylor series approximations, the extrapolation can only be 
accurate up to some fixed interval into the future, so as St is increased, N2 
must be decreased, and the speedup will decrease (of course, the accuracy of 
the full integration will then also decrease). The results shown here will change 
if Ni is changed or a different extrapolation scheme is used. A full analysis of 
projective integration for the system discussed here is beyond the scope of this 
paper (see discussions in [T2"lll4ll20f31ll32| : it is clear, however, that for problems 
with a large separation of time scales and for appropriate parameter choices, 
it will be more efficient than straight integration. Step adaptation techniques 
from traditional numerical analysis based on a posteriori error estimates can 
be modified for the adaptive selection of projective steps. 

Note that the simulations shown in Figs. [3] and H] started at t — 0, and thus 
show transient behaviour, and different initial conditions were used for the two 
simulations. It is also important to note that we chose different realisations 
of the fii at each lifting step; the results are therefore representative of the 
expected behaviour over different realizations of the random variable. Should 
we only be interested in the acceleration of computations for a particular, 
single realization, the results would be even more accurate. 



4 The 1:1 orbit 



Consider the case of 1:1 locking, i.e. solutions for which each oscillator under- 
goes one oscillation during each forcing period. The usual way to study this 
would be to "strobe" the system once each forcing cycle. Defining X p to be 
the state of the system at t = l-apjuj, where p is an integer, i.e. 

X p = X(2np) E R 2N (11) 

where X is defined in (jUJ), we could construct a map g : M. 2N — > M. 2N as 

X p+1 = g(X p ) (12) 



A 1:1 locked orbit is then a fixed point of g and its stability is determined by the 
eigenvalues of the Jacobian of g, evaluated at the fixed point. However, finding 
such a fixed point by, for example, Newton's method, is computationally very 
expensive due to the high dimensionality of the system. Also, the results we 
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obtain will only be correct for the particular realisation of the /ij, a point to 
which we return below. Instead we use the low-dimensional description of the 
system in terms of the variable Z £ IR 2 ( g+1 ). Defining Z p = Z(2np), where Z 
is defined in f llOp . we can construct a map h : R 2 ( g+1 ) — > IR 2 ( 9+1 ) as 

Z p+1 = h(Z p ) (13) 
(From now on we choose q = 1, so h : R 4 — > M 4 .) 

Note that a fixed point of h is generally not a fixed point of g\ however, we 
will see that fixed points of h do describe the overall behaviour of the system, 
and the stability follows from the eigenvalues of the Jacobian of h, evaluated 
at its fixed points. In Fig. we show the difference in x values (and in y 
values) after a time of one period, for a fixed point of h, i.e. a 1:1 locked orbit 
in the variables ao,...,b±. We can see that none of the oscillators returned 
precisely to its initial condition. The two distributions x«(0) and Xi(2n/u)) 
give the same values of and ai, even though they clearly do not completely 
coincide. Similarly for yj(0) and y^ir/uj). If the order of our approximation 
(i.e. q) was increased, the discrepancy shown in Fig. [5] would decrease and the 
fixed point of h would better approximate the fixed point of g. 

To evaluate h[Z) in practice, we lift from Z v to X p using ©-([8]), integrate ([!])- 
([2]) for one period, then restrict from X p+1 to Z p+ i using ©-([6]). Although we 
could use projective integration as described in Sec. El to integrate ([I])-©, for 
simplicity here we did not. 



5 Continuation 

We can continue fixed points of h as parameters in CO)-© are varied using 
standard pseudo-arclength continuation software [8]. In Fig. [6] we show the 1:1 
locked orbit as to is varied for a single oscillator (or equivalently, the network 
with (3 = 0, since in this case all oscillators behave identically). The left and 
right boundaries of the closed curves are saddle-node bifurcations where stable 
and unstable 1:1 locked orbits annihilate one another [7]. 

We want to analyse the case when (3 ^ 0, i.e. when the network is hetero- 
geneous. We could do this for a single realisation of the fii as above, but to 
be more general we choose a number of different realisations of the /ij and 
average over them. We do this averaging within our definition of the map h. 
Suppose r is the number of realisations we average over. For each j = 1, . . . , r 
we calculate Y p+l = h(Y p ) using the jth realisation of the /ij. (Note that Y p is 
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fixed.) We then define the averaged map h as 

kYp) = l£,YUv ( 14 ) 

' 3=1 

The results of implementing this averaging and following the 1:1 orbit are 
shown in Fig. [6] (dashed line). We can see that the effect of the heterogeneity 
is to move the range of uj values for which there is locking to lower frequencies. 
Even though the behaviour of the system was determined by following fixed 
points of h, the results agree extremely well with those found from direct 
numerical integration of the full system ([I])-© for any realisation of the /ij 
from the correct distribution. 

5. 1 Varying A 

We can follow the saddle-node bifurcations of Fig. El which mark the edges 
of the locking region, as both A and uj are varied. The results are shown in 
Fig. [7J where the resonance "tongues" for a single oscillator (solid line) and 
a network of 500 oscillators with f3 = 0.5 (dashed line) are shown. Fig. [6] is a 
horizontal slice along the top boundary of Fig. [7J 

5.2 Varying (3 

5.2.1 Breakdown of the reduced description 

It is clear that increasing (3 increases the heterogeneity of the network. To 
understand the effects of this, in Fig. [8] we plot the boundaries of the 1:1 
tongue as uj and (3 are both varied. We see the tongue boundaries move to 
lower frequencies, as expected from previous results. (Note that Fig. [HI shows 
slices through Fig. [S]at (3 = and f3 = 0.5.) When the system is unforced, (3 
and e act in opposition: if the heterogeneity (i.e. (3) is increased, the coupling 
strength (e) must be increased in order to keep the network synchronised. 
However, we consider e to be fixed. Thus for (3 large enough the forced system 
will no longer act as a "super-oscillator" in which all of the oscillators are 
synchronised with each other. Once this occurs the concept of locking between 
all oscillators and the forcing signal is no longer valid and the algorithm for 
following "coarse" (or macroscopic) saddle-node bifurcations terminates due 
to a lack of convergence within user-specified tolerances. We demonstrate this 
phenomenon in detail in Fig. 

The top two panels of Fig. show the behaviour for a typical realisation of 
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the Hi just outside the 1:1 tongue, for a high value of (3. The oscillators are 
ordered by their /i values. The behaviour of the 10 oscillators with highest ^ 
is shown in panel A. In this case, oscillators 1-492 are synchronised with each 
other, but oscillators 493-500 are not synchronised with the rest of the group. 
However, oscillators 1-492 have also lost their locking to the forcing signal, 
and this is demonstrated on panel B, where we plot X\ clS db function of time. 
This slow (apparently quasiperiodic) modulation is typical for an oscillator 
just outside a resonance tongue. (Plotting Xi for any 1 < i < 492 would give 
a similar picture.) 

Panels C and D show the behaviour just inside the tongue (note: a different 
realisation of the /i« from that in panels A and B has been used). Here, os- 
cillators 1-496 are synchronised with each other, but oscillators 497-500 are 
not synchronised with the first 496. However, now oscillators 1-496 appear to 
be still, for all practical purposes, entrained by the forcing. This is shown in 
panel D, where x\ is plotted as a function of time. The (apparently) periodic 
oscillation shown here has the same frequency as the forcing, and a plot of 
Xi for any 1 < i < 496 would be very similar. We say "apparently" periodic 
motion because once one oscillator has desynchronised from the main group, 
none of the oscillators will undergo truly periodic motion. Instead, the motion 
is expected to be quasiperiodic with at least two frequencies present, or maybe 
even weakly chaotic. 

For the results shown in Fig. [HI (with (3 = 1.2) approximately 1% of the oscil- 
lators (those with the highest values of fi) are not synchronised with the main 
cluster, either inside or outside of the tongue. However, the remaining ~ 99% 
are synchronised with each other and using the "macroscopic" approach taken 
here we can detect whether this large cluster is synchronised with the forcing 
signal or not. 

As (3 is increased, the fraction of oscillators no longer locked to the main cluster 
increases and the description of the system from the macroscopic point of 
view as a forced super-oscillator, using polynomial chaos coefficients, becomes 
increasingly flawed. This is the reason for deciding to terminate the curves in 
Fig. [HI Note that the two curves in Fig. [H] terminate at different values of f3, 
but for both curves, the saddle-node bifurcation following algorithm fails to 
converge within tolerances when approximately 1% of the oscillators become 
desynchronised from the main group. 

Note that increasing the number of Hermite polynomials, q, used in the macro- 
scopic description (thus increasing the accuracy of the low-dimensional de- 
scription) will not allow these curves to be followed to greater values of (3. It is 
the lack of synchrony within the forced network that underlies the termination 
of the curves. Of course, increasing e would allow the curves in Fig. [8] to be 
meaningfully continued to higher values of (3. 
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Note that if we were to follow a vertical path through the middle of the 
tongue shown in Fig. [8] for a particular realisation of the /ij, there would be 
many "fine-scale" bifurcations as one or more oscillators desynchronised from 
the main group. However, these are not visible in our macroscopic description 
of the system; we would need to change our macroscopic description in order 
to detect them [27J. 

Moon et al. [26|28] also considered the loss of synchrony in a heterogeneous 
network of Kuramoto oscillators. They were studying projective integration 
and showed that if one or two oscillators broke from the main cluster, projec- 
tive integration could continue, as long as the low-dimensional description was 
augmented by the phase angle (s) of the oscillators that had lost synchrony. 
We take a different approach here, regarding the unsynchronised oscillators as 
providing a perturbation to the dynamics of the synchronised group. 



5.2.2 Phase walkthrough 

For a single periodically driven oscillator, "phase walkthrough" can occur just 
outside a 1:1 resonance tongue [13J. In this phenomenon the driven oscillator 
appears to be nearly synchronised with the driving oscillator, but every so of- 
ten it undergoes either one extra or one fewer oscillation than the drive before 
returning to near synchrony. This is because the system lies in the vicinity 
of a saddle-node bifurcation of periodic orbits. This walkthrough occurs ap- 
proximately periodically, and the period scales as \u — cu*]^ 1 ^ 2 , where to* is 
the value of u> at the relevant tongue boundary [TD] . As can be seen, this slow 
oscillation can be made arbitrarily slow by adjusting uj. 

A similar phenomenon occurs in our system, but with a slight difference. 
For small (3 all of the oscillators are synchronised with one another, effectively 
acting as one oscillator, and we can observe phase walkthrough near the tongue 
boundaries with the scaling just mentioned above. However, this phenomenon 
is a result of the system spending a long time in phase space near the remains 
of the stable and unstable fixed points of h, and is thus sensitive to noise or 
other perturbations. 

Once at least one "rogue" oscillator has become desynchronised from the rest 
(as a result of increasing j3) the system can be thought of as a noisily per- 
turbed oscillator, the "oscillator" being the vast majority of oscillators that 
are synchronised with each other, and the "noise" resulting from the influence 
of the desynchronised oscillator(s) on the rest. Thus we expect that we can no 
longer make the slow oscillation arbitrarily slow just by adjusting uo. Indeed, for 
large fixed f3, near the boundaries shown in Fig. [8] there is a range of u values 
for which the slow oscillation (walkthrough) period is not well-defined, since 
perturbations from the desynchronised oscillator (s) affect the neutrally sta- 
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ble behaviour at the underlying bifurcation, resulting in apparently stochastic 
"slipping" relative to the forcing signal. It may be possible to describe these 
rare occurrences in terms of Langevin dynamics on a low-dimensional free 
energy surface [T8"f2"3"] . 



5.3 Varying (ft 

Another parameter of interest to vary is (ft. Recall that varying in a single 
unforced oscillator causes a Hopf bifurcation, leading to oscillations. The result 
of varying (ft is shown in Fig. [TU], for both a single oscillator and for a network 
with (3 = 0.3. We see that the tongue terminates at a positive value of (ft, and 
that heterogeneity moves the tongue boundary to lower values of (ft. 

To understand the cusps for low values of (ft we plot in Fig. [IT] a cross-section 
through Fig.[]JJat (ft = 0.8, for a single oscillator. The four saddle-node bifurca- 
tions are clear. (For the network, a similar plot is found, not shown.) The cusps 
involve both saddle-node bifurcations being annihilated at a codimension-two 
point. In the vicinity of these cusps, previous results on the periodically forced 
van der Pol oscillator [T7] show that there should be a curve of Hopf bifur- 
cations of the fixed point of the map starting near each cusp, which will be 
associated with the generation of quasiperiodic motion. We can follow these 
curves using standard algorithms [9], and the results for the left cusp are shown 
in Fig. [12], both for a single oscillator and for the inhomogeneous network of 
500 oscillators. (We also followed the Hopf bifurcation curve associated with 
the right cusp, not shown.) The Hopf bifurcations correspond to a complex 
conjugate pair of eigenvalues crossing out of the unit circle in the complex 
plane as u is decreased. Writing these eigenvalues at bifurcation as e ±te , we 
have 9 = at the rightmost point of the Hopf bifurcation curve (i.e. eigen- 
values of +1, +1) and 9 monotonically increases as u is decreased until 9 = ir 
(i.e. eigenvalues of —1, —1) at the leftmost point on the Hopf bifurcation curve. 

Note that while following the curve of Hopf bifurcations for the network, we 
could not use a larger value of (3 than f3 « 0.3 (i.e. we could not use a more 
heterogeneous network) because for larger values of (3 the oscillators with 
the largest values of \i would become desynchronised from the rest as the 
bifurcation was approached. The problem discussed in Sec. 15.21 regarding the 
effectiveness of the macroscopic approach would then reoccur. 

The Hopf bifurcation for a single forced oscillator is supercritical, with a stable 
2-torus being created as the Hopf bifurcation curve is crossed in the direction 
of deceasing uo [17]. The criticality of the Hopf bifurcation for the network 
seems to be the same as that for a single oscillator, and even though the curve 
in Fig. [T2] was found by averaging over 20 realisations of the /ij, it is still a 
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very good predictor of the parameter values at which quasiperiodic dynamics 
occur for any particular realisation of the /ij (not shown). 

Previous results [17] lead us to we expect (for each cusp) a curve on which 
there are orbits homoclinic to the fixed point of the map (i.e. homoclinic 
to a periodic orbit in the full system), emanating from the point where the 
curve of Hopf bifurcations and saddle-node bifurcations meet, for both a single 
oscillator and the network. We do not consider these curves further. 



6 Discussion 

In this paper we studied a finite network of heterogeneous, coupled oscilla- 
tors, all subject to the same periodic forcing. We coarse-grained the dynam- 
ics, obtaining a low- dimensional description of the system in terms of a few 
polynomial chaos coefficients. We defined a return map by sampling the low- 
dimensional system once every forcing period; by finding and following fixed 
points of this map we performed standard bifurcation analysis on the 1:1 locked 
state. By averaging over realisations of the distribution of the heterogeneity 
we have been able to obtain results valid for any particular realisation. 

We have concentrated on only the 1:1 resonance; the techniques used can be 
easily applied to any other resonances. One issue we have not discussed is 
varying N, the number of oscillators in the network. We found that for small 
values of N the number of realisations, r, of the /ij that are averaged over in 
the definition of h (eqn. (TH|) ) had to be increased in order for continuation 
algorithms to converge to a given tolerance. This makes sense since, as N is 
increased, the difference between a simulation with one particular realisation 
of the fii and that for a simulation with a different realisation, will decrease, 
and thus fewer realisations will need to be averaged over. 

Compared to the Kuramoto model results of Moon et al. [26i28j, or their 
animal flocking models [27J where each oscillator consists of a single variable 
(a phase angle) the coupled units here are representative of general ODE- 
based oscillators, capable of undergoing Hopf bifurcations. It is this extra 
feature that gives rise to results such as those in Figs. [TU] and H21 

It would be interesting to apply these ideas to networks of oscillators where the 
onset of oscillation is through different types of bifurcations, e.g. homoclinic 
bifurcation, saddle- node-on-a-circle (or SNIPER), or a saddle- node bifurcation 
of periodic orbits, or to other finite heterogeneous networks that have previ- 
ously been studied [5|30] . Another possibility is to study the periodic forcing 
of networks of coupled bursting neurons. The response of isolated bursting 
neuron models to periodic forcing has recently been studied [6TI2"2"] , as has the 
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behaviour of coupled bursters [19]. 

Regarding the problems caused by one or more oscillators desynchronising 
from the main group (see Sec. 15. 2p . one way to deal with this might be to 
expand the states of oscillators in functions other than "globally" defined 
polynomials of fi. For example, a wavelet basis of localised functions may be 
more suitable [23], particularly if the network breaks into clusters, with the 
oscillators within each cluster being synchronised. 
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Fig. 1. Angular frequency of an isolated oscillator {[3 = = e) as a function of 4>. 
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Fig. 2. Left column: function of jii (dots). Right column: t/i as Ql 

function of /ij (dots). Also included are the polynomial chaos expansions with 
q = 1 (lines). From top to bottom: t = 0,27t/oj,4it/uj. Parameters are 
A = 0.5,w = 0.85,/? = 0.1,e = l,JV = 500. 
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Fig. 3. Top: speedup (as denned in the text) as a function of N2, with N\ = 3. 
Bottom: a\ for projective integration with N2 = 1 (dots) and ai from full in- 
tegration (solid line — indistinguishable from the dots). Other parameters are 
A = 0.5, uj = 0.85, (f> = l,/3 = 0.5, e = l,q = 2. See text for details. 
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Fig. 4. Top: a\ for projective integration with N2 = 71 (dots) and a\ from full 
integration (solid line). Bottom: blowup of the small rectangle in the top panel. 
Other parameters are A = 0.5, uo = 0.85, <j> = 1,(3 = 0.5, e = l,q = 2. See text for 
details. 
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Fig. 5. Difference in x values after one period (circles), and difference in y values 
(crosses), for a fixed point of h (i.e. for an orbit that is periodic in the coefficients 
do, • • • , b\). Parameters are A = 0.5, w = 0.85, (3 = 0.5, e = 1. 
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Fig. 6. The 1:1 orbit. Solid line: x at multiples of 2tt/uj for a single oscillator. Dashed 
line: do at multiples of 2tt/lo for a network with N = 500 and f3 = 0.5, averaged 
over r = 20 realisations. Other parameters are A = 0.5, e = 1, (j> = 1. 
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Fig. 7. Boundaries of the 1:1 Arnold tongue. Solid line: one oscillator, dashed: av- 
eraging over r = 20. N = 500 and (3 = 0.5. Other parameters are e = 1, (f> = 1. 
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Fig. 8. Boundaries of the 1:1 Arnold tongue, averaging over between 20 and 60 
realisations. The curves terminate as j3 is increased because the oscillators become 
too heterogeneous to synchronise among themselves. N = 500 and A = 0.5. Other 
parameters are e = 1, (ft = 1. 
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Fig. 9. A and B: u = 0.935, j3 = 1.2 (just to the right of the right boundary shown in 
Fig. [8]). C and D: u = 0.925, (3 = 1.2 (just to the left of the right boundary). Panels 
A and C show the evolution of X491 to X500 (colour coded, white high), while panels 
B and D show x\ as a function of time. N = 500 and A = 0.5. Other parameters 
are e = 1, 6 = 1. 
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Fig. 10. Boundaries of the 1:1 orbit. Solid line: one oscillator. Dashed line: a network 
with N = 500 and A = 0.5, (3 = 0.3, r = 20, e = 1. 
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Fig. 11. x at multiples of 2tt/lo for the 1:1 orbit of a single oscillator, as u is varied. 
Solid line: stable, dashed: unstable. Parameters are A = 0.5, (ft = 0.8. This figure is 
a horizontal slice at <p = 0.8 through Fig. [TUJ 
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Fig. 12. Boundaries of the 1:1 Arnold tongue (solid and dashed lines, for a single 
oscillator and a network of N = 500 oscillators, respectively) and Hopf bifurcation 
curves (dots and circles joined by lines, for a single oscillator and a network of 500 
oscillators, respectively). Parameters are A = 0.5, /3 = 0.3, r = 20, e = 1. Compare 
with Fig. QJH 
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